sampleDists <- dist(t(assay(rld)))
plot(hclust(sampleDists))annot = dplyr::select(experimental_metadata, condition)
row.names(annot) = experimental_metadata$sample_id
rld %>%
assay() %>%
cor() %>%
pheatmap(color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
annotation = annot,
annotation_colors = list(condition = c(Control = "#BF3100", OKSM = "#E9B44C",
Senolytic = "#1B998B", Senolytic_OKSM = "#5D576B")),
cluster_rows = TRUE,
cluster_cols = T,
cellwidth = 13,
cellheight = 13)ntop = 500
rv <- rowVars(assay(rld))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca = prcomp(t(assay(rld)[select,]))
percentVar <- pca$sdev^2/sum(pca$sdev^2)
pca_data <- plotPCA(rld, intgroup = c("condition", "replicate"), returnData=TRUE)
percentVar <- round(100 * attr(pca_data, "percentVar"), digits=2)
ggplot(pca_data, aes(PC1, PC2, color=condition, shape=replicate)) + geom_point(size=3) +
scale_x_continuous(paste0("PC1: ",percentVar[1],"% variance")) +
scale_y_continuous(paste0("PC2: ",percentVar[2],"% variance")) +
coord_fixed() + theme_classic() + geom_text(data = pca_data, aes(PC1,PC2, label = name), hjust = 1.2)From the hierarchical clustering, correlation heatmap, and PCA, there seems to be some variation that we’re not accounting for. Here, we use SVA to account for the unknown variation. Since the number of variables is unknown, we call the sva function without the n.sv argument, allowing the algorithm to estimate the number of factors.
# There are batch effects for this dataset. However, running ~ batch + condition gives us an error that says the model matrix is not full rank. So we follow this guide https://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#removing-hidden-batch-effects
library(sva)
dat <- counts(dds, normalized = TRUE)
idx <- rowMeans(dat) > 1
dat <- dat[idx, ]
mod <- model.matrix(~ condition, colData(rld))
mod0 <- model.matrix(~ 1, colData(rld))
svseq <- svaseq(dat, mod, mod0)## Number of significant surrogate variables is: 3
## Iteration (out of 5 ):1 2 3 4 5
According to SVA, there are 3 significant surrogate variables.
par(mfrow = c(3, 1), mar = c(3,5,3,1))
for (i in 1:3) {
stripchart(svseq$sv[, i] ~ dds$condition, vertical = TRUE, main = paste0("SV", i))
abline(h = 0)
}After accounting for the 3 significant surrogate variables, this is what the data looks like:
sampleDists <- dist(t(assay(rld)))
plot(hclust(sampleDists))rld %>%
assay() %>%
cor() %>%
pheatmap(color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
annotation = annot,
annotation_colors = list(condition = c(Control = "#BF3100", OKSM = "#E9B44C",
Senolytic = "#1B998B", Senolytic_OKSM = "#5D576B")),
cluster_rows = TRUE,
cluster_cols = T,
cellwidth = 13,
cellheight = 13) The correlation heatmap shows that the sample AT_Sen_1 could be an outlier.
ntop = 500
rv <- rowVars(assay(rld))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca = prcomp(t(assay(rld)[select,]))
percentVar <- pca$sdev^2/sum(pca$sdev^2)
pca_data <- plotPCA(rld, intgroup = c("condition", "replicate"), returnData=TRUE)
percentVar <- round(100 * attr(pca_data, "percentVar"), digits=2)
ggplot(pca_data, aes(PC1, PC2, color=condition, shape=replicate)) + geom_point(size=3) +
scale_x_continuous(paste0("PC1: ",percentVar[1],"% variance")) +
scale_y_continuous(paste0("PC2: ",percentVar[2],"% variance")) +
coord_fixed() + theme_classic() + geom_text(data = pca_data, aes(PC1,PC2, label = name), hjust = 1.2)Taking the three plots into account, AT_Sen_1 is removed.
We re-run sva on the new dataset not including AT_Sen_1:
# There are batch effects for this dataset. However, running ~ batch + condition gives us an error that says the model matrix is not full rank. So we follow this guide https://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#removing-hidden-batch-effects
library(sva)
dat <- counts(dds, normalized = TRUE)
idx <- rowMeans(dat) > 1
dat <- dat[idx, ]
mod <- model.matrix(~ condition, colData(rld))
mod0 <- model.matrix(~ 1, colData(rld))
svseq <- svaseq(dat, mod, mod0)## Number of significant surrogate variables is: 2
## Iteration (out of 5 ):1 2 3 4 5
According to SVA, there are 2 significant surrogate variables.
par(mfrow = c(2, 1), mar = c(3,5,3,1))
for (i in 1:2) {
stripchart(svseq$sv[, i] ~ dds$condition, vertical = TRUE, main = paste0("SV", i))
abline(h = 0)
}After accounting for the 2 significant surrogate variables, this is what the data looks like:
sampleDists <- dist(t(assay(rld)))
plot(hclust(sampleDists))rld %>%
assay() %>%
cor() %>%
pheatmap(color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
annotation = annot,
annotation_colors = list(condition = c(Control = "#BF3100", OKSM = "#E9B44C",
Senolytic = "#1B998B", Senolytic_OKSM = "#5D576B")),
cluster_rows = TRUE,
cluster_cols = T,
cellwidth = 13,
cellheight = 13)ntop = 500
rv <- rowVars(assay(rld))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca = prcomp(t(assay(rld)[select,]))
percentVar <- pca$sdev^2/sum(pca$sdev^2)
pca_data <- plotPCA(rld, intgroup = c("condition", "replicate"), returnData=TRUE)
percentVar <- round(100 * attr(pca_data, "percentVar"), digits=2)
ggplot(pca_data, aes(PC1, PC2, color=condition, shape=replicate)) + geom_point(size=3) +
scale_x_continuous(paste0("PC1: ",percentVar[1],"% variance")) +
scale_y_continuous(paste0("PC2: ",percentVar[2],"% variance")) +
coord_fixed() + theme_classic() + geom_text(data = pca_data, aes(PC1,PC2, label = name), hjust = 1.2)effective_lengths = matrix(0, ncol=length(experimental_metadata$sample_id), nrow=17714)
colnames(effective_lengths)= experimental_metadata$sample_id
for( i in experimental_metadata$sample_id){
effective_lengths[,i] = read.table(paste("data/ubiquitous/", i, ".genes.results",sep=""), sep="\t", header=TRUE)$effective_length
}
row.names(effective_lengths) = read.table(paste("data/ubiquitous/", i, ".genes.results",sep=""), sep="\t", header=TRUE)$gene_id
effective_lengths = rowMeans(effective_lengths[row.names(counts(ddssva)),])
ncrpk = counts(ddssva) / (effective_lengths / 1000)
ncrpk = apply(ncrpk, c(1,2), function(x){if(is.nan(x)){0}else{x}})
ncrpk = apply(ncrpk, c(1,2), function(x){if(is.infinite(x)){0}else{x}})
ncscalingfactor = colSums(ncrpk) / 1e6
nctpm = sweep(ncrpk, 2, ncscalingfactor, "/")
nctpm.melt = melt(nctpm)
ggplot(nctpm.melt, aes(x=Var2, y=value)) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 90, colour="black", hjust = 1)) + scale_x_discrete("Sample") + scale_y_continuous("TPM")tpm.threshold = 10000
test.tpm = apply(nctpm, 1, function(x){ any(x> tpm.threshold) })
ensembl.genes[names(test.tpm[test.tpm])] %>%
as.data.frame() %>%
datatable(options = list(scrollX = TRUE))generate_de_section(dds_wald, "OKSM", "Control")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 317"
generate_de_section(dds_wald, "Senolytic", "Control")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 550"
generate_de_section(dds_wald, "Senolytic_OKSM", "Control")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 284"
generate_de_section(dds_wald, "OKSM", "Senolytic")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 460"
generate_de_section(dds_wald, "OKSM", "Senolytic_OKSM")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 170"
generate_de_section(dds_wald, "Senolytic", "Senolytic_OKSM")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 493"
dds_LRT = nbinomLRT(ddssva, reduced = ~SV1+SV2)
res = results(dds_LRT)
res$gene_biotype= ensembl.genes$gene_biotype[match(row.names(res), ensembl.genes$gene_id)]
res$external_gene_name= ensembl.genes$external_gene_name[match(row.names(res), ensembl.genes$gene_id)]
hist(res$pvalue)Number of significant genes (padj < 0.1):
sum(res$padj < 0.1, na.rm=T)## [1] 1282
# assay(x) to access the count data
assay(significant_rld) <- limma::removeBatchEffect(assay(significant_rld), covariates = cov)
sig_mat_rld = assay(significant_rld)
# The apply function swaps the rows to samples and columns to genes -- the standard is the other way around: samples in cols and genes in rows, hence the transpose function
zscores = t(apply(sig_mat_rld, 1, function(x){ (x - mean(x)) / sd(x) }))foo = as(zscores, "matrix")
bar = sapply(1:20, function(x){kmeans(foo, centers=x)$tot.withinss})
plot(bar, type="l")pam_clust <- generate_data(zscores, 7, "pam")
saveRDS(pam_clust, "output/ubiquitous/pam_clust.rds")
# pam_clust <- as.data.frame(pam_clust)
# pam_clust$Cluster <- factor(pam_clust$Cluster, levels = c(5,4,3,1,2,6))
# pam_clust <- pam_clust[order(pam_clust$Cluster),]
pheatmap(pam_clust[,1:(ncol(pam_clust)-1)],
color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
fontsize_row = 5.5,
annotation_col = annotation,
annotation_colors = anno_colours,
cluster_rows = FALSE,
cluster_cols = FALSE)pam_clust_df <- pam_clust %>%
as.data.frame() %>%
mutate(gene_name = ensembl.genes[rownames(.),]$external_gene_name) %>%
dplyr::select(gene_name, Cluster) %>%
dplyr::rename("Gene Name" = gene_name)
datatable(pam_clust_df, options = list(scrollX = TRUE), class = "white-space: nowrap")c1 <- pam_clust_df[pam_clust_df$Cluster == 1, ] %>%
dplyr::select(-Cluster)
c2 <- pam_clust_df[pam_clust_df$Cluster == 2, ] %>%
dplyr::select(-Cluster)
c3 <- pam_clust_df[pam_clust_df$Cluster == 3, ] %>%
dplyr::select(-Cluster)
c4 <- pam_clust_df[pam_clust_df$Cluster == 4, ] %>%
dplyr::select(-Cluster)
c5 <- pam_clust_df[pam_clust_df$Cluster == 5, ] %>%
dplyr::select(-Cluster)
c6 <- pam_clust_df[pam_clust_df$Cluster == 6, ] %>%
dplyr::select(-Cluster)
c7 <- pam_clust_df[pam_clust_df$Cluster == 7, ] %>%
dplyr::select(-Cluster)
data.frame(Cluster = c("Cluster 1", "Cluster 2", "Cluster 3", "Cluster 4", "Cluster 5",
"Cluster 6", "Cluster 7", "Total"),
Number_of_genes = c(nrow(c1), nrow(c2), nrow(c3), nrow(c4), nrow(c5),
nrow(c6), nrow(c7),
sum(c(nrow(c1),nrow(c2),nrow(c3),nrow(c4),nrow(c5),
nrow(c6),nrow(c7)))))## Cluster Number_of_genes
## 1 Cluster 1 300
## 2 Cluster 2 201
## 3 Cluster 3 157
## 4 Cluster 4 173
## 5 Cluster 5 147
## 6 Cluster 6 204
## 7 Cluster 7 100
## 8 Total 1282
## Checking the stability of clusters
#fviz_nbclust(zscores, kmeans, method = "silhouette")
set.seed(2)
dd = as.dist((1 - cor(t(zscores)))/2)
pam = pam(dd, 7, diss = TRUE)
sil = silhouette(pam$clustering, dd)
plot(sil, border=NA, main = "Silhouette Plot for 7 Clusters")#listEnrichrSites()
setEnrichrSite("FlyEnrichr")
dbs <- listEnrichrDbs()
to_check <- c("GO_Biological_Process_2018", "KEGG_2019")
output_dir <- "output/ubiquitous/"eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c7_go")
c1_go <- get_important_terms(eresList$GO_Biological_Process_2018)eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c7_kegg")
c1_kegg <- get_important_terms(eresList$KEGG_2019)eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c1_go")
c2_go <- get_important_terms(eresList$GO_Biological_Process_2018)eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c1_kegg")
c2_kegg <- get_important_terms(eresList$KEGG_2019)eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c2_go")
c3_go <- get_important_terms(eresList$GO_Biological_Process_2018)eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c2_kegg")
c3_kegg <- get_important_terms(eresList$KEGG_2019)eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c3_go")
c4_go <- get_important_terms(eresList$GO_Biological_Process_2018)eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c3_kegg")
c4_kegg <- get_important_terms(eresList$KEGG_2019)eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c5_go")
c5_go <- get_important_terms(eresList$GO_Biological_Process_2018)eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c5_kegg")
c5_kegg <- get_important_terms(eresList$KEGG_2019)eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c6_go")
c6_go <- get_important_terms(eresList$GO_Biological_Process_2018)eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c6_kegg")
c6_kegg <- get_important_terms(eresList$KEGG_2019)eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$GO_Biological_Process_2018, output_dir, "ubi_c4_go")
c7_go <- get_important_terms(eresList$GO_Biological_Process_2018)eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")get_important_terms_rds(eresList$KEGG_2019, output_dir, "ubi_c4_kegg")
c7_kegg <- get_important_terms(eresList$KEGG_2019)all_go = list(c1 = c1_go, c2=c2_go, c3=c3_go, c4=c4_go, c5=c5_go, c6=c6_go, c7=c7_go)
all_go <- all_go %>%
bind_rows(.id = "Cluster") %>%
group_by(Cluster) %>%
slice_min(order_by = Adjusted.P.value, n = 10)all_go %>%
mutate(Term = gsub("\\([^()]*\\)", "", Term),
Term = str_to_title(Term)) %>%
mutate(Cluster = factor(Cluster, levels = c(unique(all_go$Cluster)))) %>%
group_by(Cluster) %>%
mutate(Term = factor(Term, levels = Term)) %>%
ggplot(aes(x = Cluster, y = Term)) +
geom_point(aes(colour = Adjusted.P.value, size = Ratio)) +
ylab(NULL) +
xlab("Cluster") +
scale_color_continuous(low="red", high="blue", guide=guide_colorbar(reverse=FALSE)) +
scale_y_discrete(labels = function(x) str_wrap(x, width = 60),
limits=rev) +
#labs(fill = "Adjusted p-value") +
guides(colour = guide_colorbar(title = "Adjusted p-value", reverse = TRUE),
size = guide_legend(title = "Gene Ratio")) +
theme(axis.text.x = element_text(size=8)) +
theme_bw() +
ggtitle("GO Enrichment Analysis")all_kegg = list(c1 = c1_kegg, c2=c2_kegg, c3=c3_kegg, c4=c4_kegg, c5=c5_kegg, c6=c6_kegg, c7=c7_kegg)
all_kegg <- all_kegg %>%
bind_rows(.id = "Cluster") %>%
group_by(Cluster) %>%
slice_min(order_by = Adjusted.P.value, n = 10)all_kegg %>%
mutate(Term = gsub("\\([^()]*\\)", "", Term),
Term = str_to_title(Term)) %>%
mutate(Cluster = factor(Cluster, levels = c(unique(all_kegg$Cluster)))) %>%
group_by(Cluster) %>%
mutate(Term = factor(Term, levels = Term)) %>%
ggplot(aes(x = Cluster, y = Term)) +
geom_point(aes(colour = Adjusted.P.value, size = Ratio)) +
ylab(NULL) +
xlab("Cluster") +
scale_color_continuous(low="red", high="blue", guide=guide_colorbar(reverse=FALSE)) +
scale_y_discrete(labels = function(x) str_wrap(x, width = 60),
limits=rev) +
#labs(fill = "Adjusted p-value") +
guides(colour = guide_colorbar(title = "Adjusted p-value", reverse = TRUE),
size = guide_legend(title = "Gene Ratio")) +
theme(axis.text.x = element_text(size=8)) +
theme_bw() +
ggtitle("GO Enrichment Analysis")ame_res[[1]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[2]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[3]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[4]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[5]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[6]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[7]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))sessionInfo()## R version 4.0.5 (2021-03-31)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] sva_3.38.0 BiocParallel_1.24.1
## [3] mgcv_1.8-34 nlme_3.1-152
## [5] factoextra_1.0.7 enrichR_3.0
## [7] DT_0.17 cowplot_1.1.1
## [9] RColorBrewer_1.1-2 scales_1.1.1
## [11] reshape2_1.4.4 knitr_1.33
## [13] biomaRt_2.46.3 GenomicFeatures_1.42.3
## [15] AnnotationDbi_1.52.0 genefilter_1.72.1
## [17] DESeq2_1.30.1 SummarizedExperiment_1.20.0
## [19] Biobase_2.50.0 MatrixGenerics_1.2.1
## [21] matrixStats_0.58.0 BSgenome.Dmelanogaster.UCSC.dm6_1.4.1
## [23] BSgenome_1.58.0 rtracklayer_1.50.0
## [25] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7
## [27] Biostrings_2.58.0 XVector_0.30.0
## [29] IRanges_2.24.1 S4Vectors_0.28.1
## [31] BiocGenerics_0.36.0 forcats_0.5.1
## [33] stringr_1.4.0 dplyr_1.0.5
## [35] purrr_0.3.4 readr_1.4.0
## [37] tidyr_1.1.3 tibble_3.1.0
## [39] tidyverse_1.3.0 EnhancedVolcano_1.8.0
## [41] ggrepel_0.9.1 ggplot2_3.3.3
## [43] pheatmap_1.0.12 cluster_2.1.1
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.2.1 BiocFileCache_1.14.0
## [4] plyr_1.8.6 splines_4.0.5 crosstalk_1.1.1
## [7] digest_0.6.27 htmltools_0.5.1.1 fansi_0.4.2
## [10] magrittr_2.0.1 memoise_2.0.0 openxlsx_4.2.3
## [13] limma_3.46.0 annotate_1.68.0 modelr_0.1.8
## [16] extrafont_0.17 extrafontdb_1.0 askpass_1.1
## [19] prettyunits_1.1.1 colorspace_2.0-0 blob_1.2.1
## [22] rvest_1.0.0 rappdirs_0.3.3 haven_2.3.1
## [25] xfun_0.22 crayon_1.4.1 RCurl_1.98-1.3
## [28] jsonlite_1.7.2 survival_3.2-10 glue_1.4.2
## [31] gtable_0.3.0 zlibbioc_1.36.0 DelayedArray_0.16.3
## [34] proj4_1.0-10.1 car_3.0-10 Rttf2pt1_1.3.8
## [37] maps_3.3.0 abind_1.4-5 edgeR_3.32.1
## [40] DBI_1.1.1 rstatix_0.7.0 Rcpp_1.0.6
## [43] xtable_1.8-4 progress_1.2.2 foreign_0.8-81
## [46] bit_4.0.4 htmlwidgets_1.5.3 httr_1.4.2
## [49] ellipsis_0.3.1 farver_2.1.0 pkgconfig_2.0.3
## [52] XML_3.99-0.6 sass_0.3.1 dbplyr_2.1.1
## [55] locfit_1.5-9.4 utf8_1.2.1 labeling_0.4.2
## [58] tidyselect_1.1.0 rlang_0.4.10 munsell_0.5.0
## [61] cellranger_1.1.0 tools_4.0.5 cachem_1.0.4
## [64] cli_2.4.0 generics_0.1.0 RSQLite_2.2.6
## [67] broom_0.7.6 evaluate_0.14 fastmap_1.1.0
## [70] yaml_2.2.1 bit64_4.0.5 fs_1.5.0
## [73] zip_2.1.1 ash_1.0-15 ggrastr_0.2.3
## [76] xml2_1.3.2 compiler_4.0.5 rstudioapi_0.13
## [79] beeswarm_0.3.1 curl_4.3 ggsignif_0.6.1
## [82] reprex_2.0.0 geneplotter_1.68.0 bslib_0.2.4
## [85] stringi_1.5.3 highr_0.8 ggalt_0.4.0
## [88] lattice_0.20-41 Matrix_1.3-2 vctrs_0.3.7
## [91] pillar_1.6.0 lifecycle_1.0.0 jquerylib_0.1.3
## [94] data.table_1.14.0 bitops_1.0-6 R6_2.5.0
## [97] rio_0.5.26 KernSmooth_2.23-18 vipor_0.4.5
## [100] MASS_7.3-53.1 assertthat_0.2.1 rjson_0.2.20
## [103] openssl_1.4.3 rprojroot_2.0.2 withr_2.4.1
## [106] GenomicAlignments_1.26.0 Rsamtools_2.6.0 GenomeInfoDbData_1.2.4
## [109] hms_1.0.0 grid_4.0.5 rmarkdown_2.7
## [112] carData_3.0-4 ggpubr_0.4.0 lubridate_1.7.10
## [115] ggbeeswarm_0.6.0